#R Code for WP-18-326 

####Load Data and Libraries####
setwd("/Users/adamcasey/Documents/Toronto/Dissertation/Working Papers/WP Submission/Revision")
library(survival)
library(foreign)
library(sandwich)
library(stargazer)
library(lmtest)
library(zoo)
library(jtools)
library(ggplot2)
library(gmodels)
library(ggfortify)
library(survminer)
library(ggstance) 
library(grid)
library(dplyr)
dataorig <-read.csv("data_wp.csv", header=TRUE) #load client regime tscs data with controls
outcomes <-read.csv("diss_outcomes.csv",header=TRUE) #load regime failure data 
debruin <-read.csv("debruin.csv",header=TRUE) #load De Bruin data
seizures <-read.csv("GWF seizures.csv",header=TRUE) #Load GWF seizure group data 
vdem <-read.csv("vdemcontrols.csv",header=TRUE) #load control variables from V-Dem
jcontrols <-read.csv("joecontrols.csv",header=TRUE) #Read Data from Wright forthcoming (NAVCO, EPR)

dataorig <-rename(dataorig,
                  ally="formal_ally_patron",
                  french="fr_spons",
                  soviet="sov_spons",
                  american="us_spons",
                  russian="rus_spons",
                  british="uk_spons",
                  rev_regime ="revreg",
                  party="gwf_party",
                  personalist="gwf_personal",
                  military="gwf_military",
                  monarchy="gwf_monarch",
                  communist="sv_comm",
                  gdp_cap="e_migdppc",
                  gdp_gro="e_migdpgro",
                  mil_pers="milper",
                  mil_exp="milex",
                  oil="oil_gas_valuePOP_2014",
                  anti_gov="domestic1",
                  assass="domestic2",
                  strikes="domestic3",
                  guerrilla="domestic4",
                  crises="domestic5",
                  purges="domestic6",
                  revolution="domestic7",
                  riots="domestic8",
                  marxist_rebels="kb_ml_rebels") #rename client tscs data
outcomes <-rename(outcomes, 
                  gwf_casename ="Case",
                  sponsor="Sponsor",
                  duration="Duration",
                  rule_change="X1Rule_Change",
                  election="X2Election",
                  elec_no_inc = "X3Election_Trans",
                  uprising="X4Uprising",
                  coup="X5Coup",
                  firc="X7FIRC",
                  autogolpe="X8Autogolpe",
                  state_death ="X9State_Death",
                  american="us_spons",
                  soviet="sov_spons",
                  french="fr_spons",
                  british="uk_spons",
                  chinese="prc_spons",
                  viet="viet_spons",
                  yugo="yugo_spons",
                  russian="rus_spons",
                  egyptian="egy_spons",
                  military="mil",
                  personal="pers") #rename client cross-sec data 
#rename other dataset variables for merging 
vdem <-rename(vdem, cowcode="COWcode")
debruin <- rename(debruin, cowcode="ccode") 
seizures <-rename(seizures,
                  gwf_fail="gwf_case_fail",
                  duration="gwf_case_duration")
jcontrols <-rename(jcontrols, gwf_fail="gwf_case_fail",
                   duration="gwf_case_duration")
#Filter Datasets for relevant variables 
gwfseizures <-select(seizures, cowcode, year, gwf_country, gwf_casename, gwf_fail, gwf_fail_type,
                     duration, gwf_prior, foreignimposd, supportparty, ldr_group_military,
                     seizure_coup, seizure_rebel, seizure_foreign, seizure_uprising, seizure_election,
                     seizure_succession, seizure_family, leadermil, sectyapp_pers, paramil_pers) #select only necessary variables 
jcontrols2 <-select(jcontrols, cowcode, year, gwf_country, gwf_casename, gwf_fail, duration, gdpcap,
                    gdpcapl, oilpc, oilpcl, loggdp, logoil, mad_gdppc, mad_lgdppc) #select only necessary variables 
jcontrols3 <-filter(jcontrols2, duration!="NA") #Select only ARs
#Merge datasets 
#Merge Wright forthcoming with GWF 2018 data
jc4 <-left_join(jcontrols3, gwfseizures, row.names=TRUE)

controls <-left_join(dataorig, jc4)

controls2 <-left_join(controls, debruin)

controls3 <-left_join(controls2, vdem)

#Merge Outcomes (Cross Sectional Data) with GWF controls
outcomesGWF <-left_join(outcomes, gwfseizures)

#Filter Dataset for relevant variables
data <- select(controls3, cowcode, year, gwf_casename, startdate, enddate,
               spell, duration, gwf_fail, gwf_fail_subs, gwf_fail_type,
               gwf_fail_violent, jan1leader, duration_ldr,
               ldr_fail, ldr_fail_reg_fail, ldr_fail_reg_survive,
               spons, any_spons, satellite, spons_spell,
               duration_spons, ptcoupsuccess, ptcoupfail, regfailcoup, ldrcoup,
               american, russian, soviet, prc_spons,
               british, french, viet_spons, yugo_spons, egy_spons, comm_spons, coldwar, ally, 
               commissar, cbcount, pcount, party, personalist, military, monarchy, communist, 
               gdpcap, gdpcapl, logoil, e_civil_war, supportparty, 
               ldr_group_military, seizure_coup, seizure_rebel, seizure_election)

#Create CSV file of merged data
write.csv(data, "/Users/adamcasey/Documents/Toronto/Dissertation/Working Papers/WP Submission/Revision/finaldata.csv", row.names=TRUE)

####Functions####
#function to cluster standard errors 
vcovCluster <- function(mod, clus){ 
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(mod))!=length(clus)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(clus))
  N <- length(clus)           
  K <- mod$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K )) 
  uj  <- apply(estfun(mod), 2, function(x) tapply(x, clus, sum));
  rcse.cov <- dfc * sandwich(mod, meat = crossprod(uj)/N)
  return(coeftest(mod, rcse.cov, df.residual(mod)))
}

####Figure 1####
#Survival curves
#All Sponsorship 
surv1 <- survfit(Surv(duration,duration+1,gwf_fail) ~ spons, data=data)
#US sponsorship
surv2 <-survfit(Surv(duration,duration+1,gwf_fail) ~ american, data=data)
#Soviet sponsorship 
surv3 <-survfit(Surv(duration,duration+1,gwf_fail) ~ soviet, data=data)
#French sponsorship 
surv4 <-survfit(Surv(duration,duration+1,gwf_fail) ~ french, data=data)

#plot all survival curves together 
splots <-list()
splots [[1]] <- ggsurvplot(surv1, data=data,
                           palette =c("black","darkgray"),
                           linetype=c("solid","longdash"),
                           risk.table=FALSE,
                           conf.int = TRUE, 
                           xlim=c(0,110),
                           break.time.by =10, 
                           ggtheme=theme_classic(base_size=11),
                           risk.table.y.text=FALSE,
                           censor.size =1,
                           title="Survival Curves",
                           subtitle="Autocratic Regimes, 1946-2010",
                           legend.title="",
                           legend.labs=c("Other Regimes","Client Regimes"))
splots[[2]] <- ggsurvplot(surv2, data=data,
                          palette =c("black","darkgray"),
                          linetype=c("solid","longdash"),
                          risk.table=FALSE,
                          conf.int = TRUE, 
                          xlim=c(0,110),
                          break.time.by =10, 
                          ggtheme=theme_classic(base_size=11),
                          risk.table.y.text=FALSE,
                          censor.size =1,
                          legend.title="",
                          legend.labs=c("Other Regimes","U.S. Client Regimes"))
splots[[3]] <- ggsurvplot(surv3, data=data,
                          palette =c("black","darkgray"),
                          linetype=c("solid","longdash"),
                          risk.table=FALSE,
                          conf.int = TRUE, 
                          xlim=c(0,110),
                          break.time.by =10, 
                          ggtheme=theme_classic(base_size=11),
                          risk.table.y.text=FALSE,
                          censor.size =1,
                          title="",
                          subtitle="",
                          legend.title="",
                          legend.labs=c("Other Regimes","Soviet Client Regimes"))
splots[[4]] <- ggsurvplot(surv4, data=data,
                          palette =c("black","darkgray"),
                          linetype=c("solid","longdash"),
                          risk.table=FALSE,
                          conf.int = TRUE, 
                          xlim=c(0,110),
                          break.time.by =10, 
                          ggtheme=theme_classic(base_size=11),
                          risk.table.y.text=FALSE,
                          censor.size =1, 
                          legend.title="",
                          legend.labs=c("Other Regimes","French Client Regimes"))
arrange_ggsurvplots(splots, print=TRUE,
                    ncol=2, nrow=2)

##### Table 1 #####
#Regime survival cox 
#Model 1
coxreg1 <-coxph(Surv(duration,duration+1, gwf_fail) ~ american + soviet + french + british,
                method="efron",data=data)
#Model 2 
coxreg2 <-coxph(Surv(duration,duration+1,gwf_fail) ~ american + soviet + french + british + 
                    gdpcapl + logoil + e_civil_war, method="efron", data=data)
#Model 3
coxreg3 <-coxph(Surv(duration,duration+1,gwf_fail) ~ american + soviet + french + british + 
                    supportparty + seizure_coup, method="efron", data=data)
#Model 4
coxreg4 <-coxph(Surv(duration,duration+1,gwf_fail) ~ american + soviet + french + british + 
                    seizure_coup + seizure_rebel + seizure_election, method="efron", data=data)


##### Figure 2 #####
#Plot logit regime ending military coups
#Model 1
dataregcoup1 <-na.omit(data[,c('cowcode','regfailcoup','spons')])
dataregcoup1 <-rename(dataregcoup1, sponsorship=spons)
logregcoup1 <- glm(regfailcoup ~ sponsorship, 
                   data=dataregcoup1, family="binomial")
seregcoup1 <- vcovCluster(logregcoup1, dataregcoup1$cowcode)

#Model 2
dataregcoup2 <-na.omit(data[,c('cowcode','regfailcoup','american')])
logregcoup2 <-glm(regfailcoup ~ american,
                  data=dataregcoup2, family="binomial")
seregcoup2 <-vcovCluster(logregcoup2, dataregcoup2$cowcode)

#Model 3
dataregcoup3 <- na.omit(data[,c('cowcode','regfailcoup',
                                'french')])
logregcoup3 <-glm(regfailcoup ~ french,
                  data=dataregcoup3, family="binomial")
seregcoup3 <-vcovCluster(logregcoup3, dataregcoup3$cowcode)

#Model 4
dataregcoup4 <-na.omit(data[,c('cowcode','regfailcoup','british')])
logregcoup4 <-glm(regfailcoup ~ british,
                  data=dataregcoup4, family="binomial")
seregcoup4 <-vcovCluster(logregcoup4, dataregcoup4$cowcode)

#Generate plot
brc1 <- plot_coefs(logregcoup1, logregcoup2,
                   logregcoup3,
                   logregcoup4,
                   ci_level=0.95,
                   inner_ci_level=0.90,
                   colors=c("grey","black","darkgrey","lightgrey"),
                   legend.title="Model", 
                   title="Logistic Models",
                   subtitle="Regime-Ending Military Coups", 
                   scale=TRUE, robust=TRUE)
ggpar(brc1, title="Logistic Models",
      subtitle="Regime-Ending Military Coups")



##### Figure 3 #####
#Plot logit U.S. sponsorship regime ending miltiary coups 

#Model 1
dataregcoup3.1 <-na.omit(data[,c('cowcode','regfailcoup','american')])
logregcoup3.1 <-glm(regfailcoup ~ american,
                    data=dataregcoup3.1, family="binomial")
seregcoup3.1 <-vcovCluster(logregcoup3.1, dataregcoup3.1$cowcode)

#Model 2
dataregcoup3.2 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'supportparty')])
logregcoup3.2 <-glm(regfailcoup ~ supportparty,
                    data=dataregcoup3.2, family="binomial")
seregcoup3.2 <-vcovCluster(logregcoup3.2, dataregcoup3.2$cowcode)

#Model 3
dataregcoup3.3 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'seizure_coup')])
logregcoup3.3 <-glm(regfailcoup ~ seizure_coup,
                    data=dataregcoup3.3, family="binomial")
seregcoup3.3 <-vcovCluster(logregcoup3.3, dataregcoup3.3$cowcode)

#Model 4
dataregcoup3.4 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'ldr_group_military')])
logregcoup3.4 <-glm(regfailcoup ~ ldr_group_military,
                    data=dataregcoup3.4, family="binomial")
seregcoup3.4 <-vcovCluster(logregcoup3.4, dataregcoup3.4$cowcode)

#Model 5
dataregcoup3.5 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'american','supportparty')])
logregcoup3.5 <-glm(regfailcoup ~ american + supportparty,
                    data=dataregcoup3.5, family="binomial")
seregcoup3.5 <-vcovCluster(logregcoup3.5, dataregcoup3.5$cowcode)

#Model 6
dataregcoup3.6 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'american','seizure_coup')])
logregcoup3.6 <-glm(regfailcoup ~ american + seizure_coup,
                    data=dataregcoup3.6, family="binomial")
seregcoup3.6 <-vcovCluster(logregcoup3.6, dataregcoup3.6$cowcode)

#Model 7
dataregcoup3.7 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'american','ldr_group_military')])
logregcoup3.7 <-glm(regfailcoup ~ american + ldr_group_military,
                    data=dataregcoup3.7, family="binomial")
seregcoup3.7 <-vcovCluster(logregcoup3.7, dataregcoup3.7$cowcode)

#Generate coefficient plot
mrc1 <- plot_coefs(logregcoup3.1, logregcoup3.2,
                   logregcoup3.3, 
                   logregcoup3.4, logregcoup3.5,
                   logregcoup3.6, logregcoup3.7,
                   ci_level=0.95,
                   inner_ci_level=0.90,
                   colors=c("black","black","black","black", "grey","grey","grey"),
                   legend.title="Model", 
                   title="Logistic Models",
                   subtitle="Regime-Ending Military Coups", 
                   scale=TRUE, robust=TRUE)
ggpar(mrc1, title="Logistic Models",
      subtitle="Regime-Ending Military Coups")


#### Figure 4 #####
#Plot logit French sponsorship regime ending military coups 
#Model 1
dataregcoup4.1 <-na.omit(data[,c('cowcode','regfailcoup','french')])
logregcoup4.1 <-glm(regfailcoup ~ french,
                    data=dataregcoup4.1, family="binomial")
seregcoup4.1 <-vcovCluster(logregcoup4.1, dataregcoup4.1$cowcode)

#Model 2
dataregcoup4.2 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'supportparty')])
logregcoup4.2 <-glm(regfailcoup ~ supportparty,
                    data=dataregcoup4.2, family="binomial")
seregcoup4.2 <-vcovCluster(logregcoup4.2, dataregcoup4.2$cowcode)

#Model 3
dataregcoup4.3 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'seizure_coup')])
logregcoup4.3 <-glm(regfailcoup ~ seizure_coup,
                    data=dataregcoup4.3, family="binomial")
seregcoup4.3 <-vcovCluster(logregcoup4.3, dataregcoup4.3$cowcode)

#Model 4
dataregcoup4.4 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'ldr_group_military')])
logregcoup4.4 <-glm(regfailcoup ~ ldr_group_military,
                    data=dataregcoup4.4, family="binomial")
seregcoup4.4 <-vcovCluster(logregcoup4.4, dataregcoup4.4$cowcode)

#Model 5
dataregcoup4.5 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'french','supportparty')])
logregcoup4.5 <-glm(regfailcoup ~ french + supportparty,
                    data=dataregcoup4.5, family="binomial")
seregcoup4.5 <-vcovCluster(logregcoup4.5, dataregcoup4.5$cowcode)

#Model 6
dataregcoup4.6 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'french','seizure_coup')])
logregcoup4.6 <-glm(regfailcoup ~ french + seizure_coup,
                    data=dataregcoup4.6, family="binomial")
seregcoup4.6 <-vcovCluster(logregcoup4.6, dataregcoup4.6$cowcode)

#Model 7
dataregcoup4.7 <-na.omit(data[,c('cowcode','regfailcoup',
                                 'french','ldr_group_military')])
logregcoup4.7 <-glm(regfailcoup ~ french + ldr_group_military,
                    data=dataregcoup4.7, family="binomial")
seregcoup4.7 <-vcovCluster(logregcoup4.7, dataregcoup4.7$cowcode)

#Generate plot
mrc2 <- plot_coefs(logregcoup4.1, logregcoup4.2,
                   logregcoup4.3, 
                   logregcoup4.4, logregcoup4.5,
                   logregcoup4.6, logregcoup4.7,
                   ci_level=0.95,
                   inner_ci_level=0.90,
                   colors=c("black","black","black","black", "grey","grey","grey"),
                   legend.title="Model", 
                   title="Logistic Models",
                   subtitle="Regime-Ending Military Coups", 
                   scale=TRUE, robust=TRUE)
ggpar(mrc2, title="Logistic Models",
      subtitle="Regime-Ending Military Coups")


####Table 3####
#Model 1
datacpA2.1 <-na.omit(data[,c('cowcode','spons','cbcount')])
lmcpA2.1 <- lm(cbcount ~ spons, data=datacpA2.1)
selmcpA2.1 <-vcovCluster(lmcpA2.1, datacpA2.1$cowcode)
#Model 2
datacpA2.2 <-na.omit(data[,c('cowcode','american',
                             'french','soviet','british','cbcount')])
lmcpA2.2 <- lm(cbcount ~ american + soviet + french + british, data=datacpA2.2)
selmcpA2.2 <-vcovCluster(lmcpA2.2, datacpA2.2$cowcode)
#Model 3
datacpA2.3 <-na.omit(data[,c('cowcode','american',
                             'french','soviet','british','cbcount',
                             'ldr_group_military','supportparty')])
lmcpA2.3 <- lm(cbcount ~ american + soviet + french + british +
                 ldr_group_military + supportparty, data=datacpA2.3)
selmcpA2.3 <-vcovCluster(lmcpA2.3, datacpA2.3$cowcode)
#Model 4
datacpA2.4 <-na.omit(data[,c('cowcode','american',
                             'french','soviet','british','cbcount',
                             'seizure_coup','seizure_rebel','seizure_election')])
lmcpA2.4 <- lm(cbcount ~ american + soviet + french + british +
                 seizure_coup + seizure_rebel + seizure_election, data=datacpA2.4)
selmcpA2.4 <-vcovCluster(lmcpA2.4, datacpA2.4$cowcode)

#### Figure 6 #####
#Plot logit coup proofing and coups 
#Model 1
dataregcoup6.1 <-na.omit(data[,c('cowcode','regfailcoup','spons')])
dataregcoup6.1 <-rename(dataregcoup6.1, sponsorship=spons)
logregcoup6.1 <- glm(regfailcoup ~ sponsorship, 
                     data=dataregcoup6.1, family="binomial")
seregcoup6.1 <- vcovCluster(logregcoup6.1, dataregcoup6.1$cowcode)
#Model 2
dataregcoup6.2 <-na.omit(data[,c('cowcode','regfailcoup','cbcount')])
logregcoup6.2 <- glm(regfailcoup ~ cbcount,
                      data=dataregcoup6.2, family="binomial")
seregcoup6.2 <-vcovCluster(logregcoup6.2, dataregcoup6.2$cowcode)
#Model 3
dataregcoup6.3 <-na.omit(data[,c('cowcode','regfailcoup','commissar')])
logregcoup6.3 <-glm(regfailcoup ~ commissar,
                    data=dataregcoup6.3, family="binomial")
seregcoup6.3 <- vcovCluster(logregcoup6.3, dataregcoup6.3$cowcode)
#Model 4
dataregcoup6.4 <-na.omit(data[,c('cowcode','regfailcoup','spons','cbcount')])
dataregcoup6.4 <-rename(dataregcoup6.4, sponsorship=spons)
logregcoup6.4 <- glm(regfailcoup ~ sponsorship + cbcount,
                     data=dataregcoup6.4, family="binomial")
seregcoup6.4 <- vcovCluster(logregcoup6.4, dataregcoup6.4$cowcode)
#Model 5
dataregcoup6.5 <-na.omit(data[,c('cowcode','regfailcoup','spons','commissar')])
dataregcoup6.5 <-rename(dataregcoup6.5, sponsorship=spons)
logregcoup6.5 <- glm(regfailcoup ~ sponsorship + commissar,
                     data=dataregcoup6.5, family="binomial")
seregcoup6.5 <- vcovCluster(logregcoup6.5, dataregcoup6.5$cowcode)

#Generate plot
mcp1 <-plot_coefs(logregcoup6.1, logregcoup6.2, 
                  logregcoup6.3, logregcoup6.4,
                  logregcoup6.5,
                  colors=c("black","black","black","grey","grey"), 
                  ci_level=0.95,
                  inner_ci_level=0.90,
                  title="Logistic Models",
                  subtitle="Regime-Ending Military Coups", 
                  scale=TRUE, robust=TRUE)
ggpar(mcp1, title="Logistic Models",
      subtitle="Regime-Ending Military Coups")
